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1.  Introduction 


Concepts  from  singular  perturbations  have  recently  suggested  very 
useful  new  techniques  for  reduced  order  modeling,  stability  analysis,  and 
synthesis  of  optimal  controls  (cf.,  e.g.,  Kokotovic  et  al.  (1976)  and 
O'Malley  (1978)).  Their  primary  advantage  is  achieved  by  approximating 
solutions  of  systems  of  the  form 


(1) 


[ . 

u = f (u,  v,  t,  e) 

< 

ev  = g(u,v,t,e) 


with  solutions  of  the  lower  order  ("reduced"  or  "outer") 


(2) 


t 

U = f(U,V,t,0) 

< 

0 = g(U,V,t,0) 


system 


m/tcEo 


While  Sectioi 
Scctiai 


obtained  by  setting  the  small  positive  parameter  e equal  to  zero.  In 
engineering  contexts,  analogous  procedures  have  commonly  (and  necessarily) 
been  used  to  lower  the  (often  prohibitively  high)  dimensionality  of  complex 
models.  These  are  often  explained  in  terms  of  neglecting  fast  as  opposed 
to  slow  dynamics  (cf.  Davison  (1966),  van  Ness  (1977),  and  Skira  and  De  Hoff 
(1977)).  Despite  remarkable  success,  these  schemes  sometimes  involve  nonsen- 
sical use  of  asymptotic  analysis,  valid  as  e ■*  0,  for  e = 1 (cf.  Calise 
(1976)).  Indeed,  such  success,  so  perilously  based,  requires  more  careful 
analysis.  In  the  numerical  methods  literature  such  "singularly  perturbed” 


□ r 
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systems  of  differential  equations  are  called  stiff  (cf.  Willoughby  (1974)). 
Singular  perturbations  theory  has  recently  contributed  to  the  analysis  of 
such  stiff  equations  (cf.,  e.g.,  Miranker  (1975),  Hemker  (1977),  Flaherty 
and  O'Malley  (1977),  and  Kreiss  (1978)).  In  physical  contexts,  however,  one 
is  seldom  presented  with  systems  in  the  form  (1).  The  fundamental  question 
usually  remains,  viz:  How  does  one  numerically  identify  the  small  param- 
eter (s)  involved?  This  must  be  answered  before  one  can  effectively 
utilize  the  singular  perturbations  machinery,  either  directly  or  as  an 
underpinning  structure  for  a numerical  approximations  procedure. 

As  the  simplest  example  of  such  problems,  let  us  consider  the  linear 
constant  coefficient  system 

(3)  x = Ax 

on  a finite  interval,  say  0 £ t £ 1,  when  the  spectrum  A(A)  of  A has 
two  (or  more)  time  scales  (i.e.,  the  eigenvalues  of  A cluster  into  two 

or  more  sets  which  are  widely  separated  in  magnitude).  Specifically,  we 

shall  take  x to  be  an  n-vector  with  the  spectrum  of  A being  the 
disjoint  sum 

(4)  A (A)  = S U F 

where 

S = (s^,,  S2»  . } and  | s_J  £ lsi+ll  f°r  eac^ 

F = {f  , f , ....  f } with  | f . | < If.,,  | for  each  j, 

12  n_  j — j+1 
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i 


+ n2  = n with  0 < n^  < n,  and 


Sn  « fl 

nl  1 


Thus,  a small  parameter  appropriate  for  our  analysis  of  (3)  has  been 
identified  as 


e = |sn  /f,l  « 1. 

nl  1 


More  refined  partitioning  of  1(A)  (corresponding  to  several  time  scales) 
will  often  be  useful  as  well,  providing  several  small  parameters  of  de- 
creasing size.  Our  general  approach  may  even  be  useful  when  e is  not  so 
very  small,  as  we've  demonstrated  numerically. 

We  shall  also  assume  that 


I f j I >>  I sn  l»  j=l»  •••»  • 


This  will  allow  us  to  approximate  the  solutions  of  (3)  by  solutions  of  a 
lower  order  (n^-dimensional)  system,  except  in  narrow  endpoint  boundary 
layers.  We'd  have  to  consider  highly  oscillatory  solutions  and  use  aver- 
aging methods  if  we  omitted  assumption  (7)  and  allowed  some  eigenvalues 
fj  to  be  purely  imaginary  (or  nearly  so)  (cf.  Hoppensteadt  and  Miranker 
(1976)  and  Chow  (1977),  however). 


2.  Preliminary  Mathematical  Analysis 

The  solutions  of  (3)  can,  of  course,  be  obtained  by  constructing  the 
At  At 

matrix  exponential  e , or  e K for  any  nonsingular  matrix  K.  As 


-5- 


(12) 


A(B,)  - F - (f.,  ....  f }. 
2 1 


We  note  that  the  transformation  (9)  is  suggested  by  that  commonly  used 
for  time-varying  systems  (cf.,  e.g.,  Harris  (1973)).  Setting  y 1 
we  have  the  separate  systems 


(13) 


Vl 


and 


(14) 


Now  the  singular  perturbations  nature  of  (10),  or  equivalently  (13)- (14), 
would  be  made  even  more  obvious  if  we  had  multiplied  (14)  through  by  the 
small  parameter  II  II  / II  8 and  compared  the  result  with  (1). 

We  shall  obtain  the  block  diagonal  system  (10)  in  two  steps.  If  we 


first  take  K = 0 in  (9),  we  find  that 
upper  block  triangular  system 


Ly2  J 


satisfies  the  decoupled 


(15) 


provided  the  (non-square)  decoupling  matrix  L satisfies  the  algebraic 
Riccati  equation 


(16)  ^11  " A22L  ‘ M12L  + A21  = ° 

in  addition  to  the  eigenvalue  placement  requirements  (11)- (12).  (We  note 


6- 


that  (16)  is  linear  when  A^2  = 0.  Then,  however,  A is  already  in 
block-triangular  form,  although  not  generally  decoupled  into  slow  and  fast 
modes.)  We  next  obtain  the  diagonal  system  (10)  by  requiring  K to  sat- 
isfy the  linear  equation 

(17)  KB2  - B2K  + A12  = 0. 

This  Liapunov  equation  will  have  a unique  solution  since  B^  and  B2 
have  no  eigenvalues  in  common  (cf.,  e.g.,  Bellman  (1970)).  An  attractive 
feature  of  the  nonsingular  transformation  (9)  is  that  its  inverse  is 
simply  given  by 


(18) 


y 


I + KL 
nl 

L 


so  we  don't  need  to  numerically  perform  a matrix  inversion  to  solve  (9) 
for  y. 

Presuming  we  can  solve  for  L and  K,  we  can  treat  the  system  (3) 

as  the  transformed  system  (10)  and  then  separately  integrate  (13)  for 

the  slowly- varying  vector  y^  and  (14)  for  the  rapidly-changing  vector 

y„.  If  m of  the  (large)  eigenvalues  f.  of  B have  negative  real 
1.  J A 

parts,  there  will  be  an  m^-dimensional  manifold  of  solutions  to  (14) 
which  decay  rapidly  to  zero  away  from  t = 0 and  an  m2  = n2  ” mi  dimen- 
sional manifold  of  rapidly  growing  solutions  there.  Alternatively, 

B2(t-1) 

e has  an  m2  dimensional  column  space  of  vectors  which  decay 

rapidly  to  zero  within  [0,1),  away  from  t = 1.  Let  us  use  the  (not 
necessarily  Jordan)  decomposition 


* 

6 

r 

V 


i- ' 


\ f ■ 
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(19) 


B2  - ND1I 


-1 


where  D has  the  block  diagonal  form 


°1  ° 


0 D, 


with  and  -D^  being  and  x dimensional  stable 

matrices,  and 


N = N2] 


where  is  n2  * m^  dimensional.  Then 


V V 

e N1  = Nie 


spans  all  solutions  of  (14)  which  decay  away  from  t = 0,  while 


B (t-1)  D (t-1) 

e N2  = N2e  Z 


spans  all  solutions  of  (14)  which  decay  away  from  t = 1.  For  the  matrix 


-8- 


will  be  a bounded  fundamental  matrix  for  the  system  (10).  Our  use  of  the 

V 

e notation  is  harmless,  since  any  standard  numerical  code  or  explicit 

formulas  for  the  matrix  exponential  will  serve  well  to  produce  this  slowly- 

varying  solution  to  the  matrix  version  of  (13) . Unless  B or  -B^  is 

stable,  however,  e is  very  difficult  to  compute,  since  it  contains 

both  rapidly  growing  and  rapidly  decaying  modes  (though  no  slow  modes) . 

Less  difficulty  is  involved  in  separately  computing  its  fast  initial  and 

B t B (t-1) 

terminal  transients  e and  e N^.  These  correspond  to  inte- 

grating (14)  on  the  appropriate  m^  or  m ^ dimensional  manifold.  If  we 
wished  to  obtain  a bounded  solution  to  the  initial  value  problem  for  (3), 
to  obtain  a bounded  solution  on  [0,1]  as  e ->  0,  we'd  need  m^  = n^ 
or  we'd  have  to  restrict  the  initial  vector  x(0)  to  a lower  dimensional 
manifold.  We  might  even  consider  such  problems  on  t >_  0 if  A were 
stable.  Two-point  problems,  of  course,  need  not  have  a solution,  but 
we've  at  least  outlined  a framework  in  which  we  might  numerically  consider 
the  question.  These  considerations,  of  course,  parallel  those  always 
encountered  in  using  fundamental  matrices  (cf.,  e.g.,  Coppel  (1965)  and 
Cole  (1968)). 

Bt 

It  is  convenient  to  decompose  the  fundamental  matrix  e C for  (10) 
as 

(21)  V - eBtC  - 0ts  Yf()  Yfl) 
with  n^  slow  modes 

(22) 


-9- 


fast-decaying  initial  modes 


/ Y 


(23) 


fO 


fOl 


\ Yf 02 


I 0 

V 

e N, 


and  m^  fast-decaying  terminal  modes 


(24) 


fl 


1 Yfll 

i 0 

= 

= 

B (t-1) 

Yf  12  / 

l • N2 

The  decomposition  (21)  implies  a corresponding  decomposition  for  a bounded 


At 


fundamental  matrix  e C for  our  original  system  (3),  i.e., 


(25) 


X - eAcc~  i (Xs  Xf0  xfl) 


with  n^  slow  modes 


(26) 


X = 
s 


1 Xsl  1 



" 1 

Y = 

I 1 

ni 

Xs2  1 

1 + LK 

n0  / 

s 

-L 

V 


and  n^  fast  modes 


(27) 


fO 


/ xfoi 

-K  \ 

l Xf02  1 

V + LK' 

V 

e N, 


and 
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(28) 


X 

I xfll  1 

| 

1 “ ) 

fl  1 

1 Xfl2  1 

l I + LK 
\ n2  I 

B (t-1) 
e N„ 


(If  m2  = 0,  e.g.,  is  omitted  in  (25).)  Within  (0,1),  then,  all 

solutions  to  (3)  are  nicely  approximated  in  the  column  space  of  the 

reduced  order  n * n1  matrix  X (t) . In  the  initial  "boundary  layer," 

■A-  s 

i.e.,  in  a narrow  interval  near  t - 0,  X is  nearly  constant,  so  solu- 

s 

tions  are  well  approximated  by  the  span  of  the  n * (n^  + m^)  matrix 

(X  (0)  Xcri) , while  (X  (1)  XC1)  is  appropriate  near  t = 1.  If  we're  not 

concerned  with  the  fast  endpoint  transients,  we'd  simply  use  the  reduced 

(or  "outer")  approximation  X (t)  everywhere.  Since  X (0)  has  rank 

s s 

n^  < n,  while  x may  be  determined  by  n linearly  independent  endcondi- 
tions,  we  realize  that  this  approximation  is  usually  not  correct  at  t = 0 
and  1.  The  approximation  of  x by  a vector  multiple  of  Xg,  however, 
improves  as  the  ratio  e decreases  (i.e.,  the  stiffer  our  system  (3),  the 
more  negligible  are  the  fast  modes  X^q  and  X^  within  (0,1)).  We  note 
that  the  separai ion  of  growing  and  decaying  modes  is  common  in  time-varying 
stability  theory  and  is  formalized  in  the  concept  of  an  exponential  dichot- 
omy (cf.  Coppel  (1978)).  It  has  previously  been  used  by  Ferguson  (1975) 
for  the  numerical  solution  of  singular  perturbation  problems.  Finally,  we 
emphasize  the  numerical  importance  of  decoupling  the  slow  part  and  the  fast 
initial  and  terminal  transients.  Large  integration  steps  can  be  used 
throughout  to  calculate  Xg , while  an  accurate  integration  of  X^  and 
X^  (with  much  smaller  stepsizes)  can  be  obtained,  if  needed,  on  the 
appropriate  short  t- intervals.  We  intend  to  primarily  emphasize  the 
numerical  aspects  involved  with  computing  the  reduced  order  model  Xg. 

Readers  should  note  that  this  primarily  involves  the  decoupling  matrix  L, 


-11- 


alt  hough  the  matrix  K is  needed  to  identify  the  endconditions  appropriate 

for  X . 
s 


3.  The  Decoupling  Matrices  L and  K 

In  or’er  to  utilize  the  structure  we've  developed,  we  must  be  able  to 
find  an  n ^ x n^  matrix  solution  L of  the  algebraic  Riccati  equation 
(16)  satisfying  the  eigenvalue  placement  conditions  (11)- (12)  as  well.  We 
note  that  the  well-studied  case  when  A is  a Hamiltonian  matrix  is  very 
important  in  applications  (then  = -A^j  A12  = A12’  A21  = A21 

the  prime  denoting  transposition),  and  the  solution  of  (16)  is  symmetric). 
Under  assumptions  of  stabilizability  and  detectability,  (16)  has  a unique 
positive  semi-definite  solution  (cf.  Kucera  (1973)).  An  analogous  approach 
is  required  for  our  problem.  Specifically,  let  us  suppose  that  A has  the 
decomposition 


(29) 


A = MJM 


where 


(30) 


with  the  blocks  J and  J. 

s f 


satisfying 


(31)  A(J  ) = A(B.)  = S,  A(J.)  = A (B  ) = F, 

si  I l 


and  the  matrix  M having  the  compatible  partitioning 


11 

M12 

21 

M22 

(32) 


M = 
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We  could  take  M to  be  a modal  matrix  so  that  J is  in  Jordan  form.  Then 
the  columns  of  M would  be  right  eigenvectors  and  generalized  right  eigen- 
vectors of  A corresponding  to  the  diagonal  entries  of  J,  while  the  rows 
of  M-1  would  be  corresponding  (sometimes  generalized)  left  eigenvectors 
of  A (cf.,  e.g.,  Hirsch  and  Smale  (1974)).  Tt's  not  necessary,  however, 
to  be  so  explicit.  All  that  is  needed  is  the  slow-fast  block  decomposition 
of  J provided  by  (30)- (31).  Indeed,  the  Jordan  form  is  not  particularly 

well  suited  to  numerical  computation  (cf.  Golub  and  Wilkinson  (1976)). 

At  Jt 

Since  e M = Me  is  a fundamental  matrix  for  (3),  it  follows 
that  the  n^-dimensional  space  of  slowly-varying  solutions  of  (3)  coincides 
with  the  column  span  of 


(33) 


M 


J t 
s 
e 

0 

M11 

0 

0 

_M2!_ 

J t 
s 


V 

e being  fast-varying.  Comparing  with  the  slowly-varying  matrix  X ci 


1 \ 1 

( M11  \ 

(26),  it  follows  that 

and  | 

\ -L  I 

l M21  / 

Therefore,  if  we  have  a solution  L to  (16), 
invertible  and  we  must  have 


must  span  the  same  space. 

(11),  and  (12),  Mn  will  be 


(34) 


L = 


-M  M . 
2l"ll 


Conversely,  if  LM.^  = -M21  ^or  M11  invertible,  multiplying  out 


I 

0 

I 

0 

nl 

A 

nl 

L 

I 

-L 

I 

n2_ 

n2_ 
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and  using  the  decomposition  (29)  for  A produces  the  lower  right  entry 
LA^  - A^L  - LA^L  + A21  = ® (after  considerable  manipulation).  Further, 
if 


___ 

__ 

Q11 

CM 

rH 

O' 

(35) 

M_1  = 

_Q21 

q22_ 

we’ll  have  Q22  = (M22  - M2;]M~Jm  ) 1 and  (34)  will  hold  if  and  only  if 


Kv  . , 


(36) 


L ■ ^2Q21' 


The  formulas  (34)  and  (36)  show  that  the  matrix  L we  are  seeking  is 
unique,  although  M is  far  from  unique  and  the  quadratic  equation  (16) 
has  many  other  solutions.  Wu  and  Narasimhamurthi  (1976),  Narasimhamurthi 
and  Wu  (1977),  and  Anderson  (1978)  contain  alternative  and  more  detailed 
treatments  of  these  results. 

In  order  to  compute  L,  we  first  employ  a standard  eigenanalysis 
library  routine  to  obtain  approximate  eigenvalues  of  the  matrix  A.  This 
allows  us  to  select  n , the  number  of  moderate  eigenvalues,  and  thereby 
determines  the  number  n2  = n - n^  of  eigenvalues  with  large  real  parts. 

We  note  that  multiple  and  complex  conjugate  eigenvalues  are  naturally 
grouped  together  in  either  S or  F. 

If  n^  £ n/2,  we  obtain  n^  approximate  right  (or  generalized 
right)  eigenvectors  corresponding  to  the  n^  moderate  eigenvalues,  thereby 


forming  an  approximation 


rM°  n 

M ~ 1 

ii 

11 

to  the  submatrix 

M° 

_21_ 

_M2!_ 

of  (32).  (If 


the  resulting  matrix  has  complex  entries,  we  eliminate  them  through 
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elementary  column  operations.)  Then  we  solve  the  resulting  real  system 


(37) 


LoMli 


approximating  (34),  providing  M°^  is  nonsingular.  An  efficient  method 
of  obtaining  an  approximate  solution  to  (37)  is  the  "LU"  factoriza- 

tion technique  (cf.,  e.g.,  Stewart  (1973)).  If  n^  > n/2,  it  would  be 
more  efficient  to  seek  n^  < n/2  approximate  left  (or  generalized  left) 
eigenvectors  of  A and  instead  obtain  an  approximate  solution  to  (36). 

In  either  case,  we  don't  need  all  the  eigenvectors  of  A.  For  all  our 
computations  to  date,  we  have  used  a modal  matrix  for  M,  but  we  expect 
to  not  be  so  restrictive  in  later  testing. 

Thus  far,  we  have  used  x^  and  x^  as  any  n^  and  n^  dimen- 


sional subvectors  of  x,  requiring  only  that  the  matrix 


M 

M 


11 

21 


of 


slow  modes  (necessarily  of  rank  n^)  have  nonsingular.  This  may 

demand  a reordering  of  the  components  of  x.  If  we  used  an  ordering  for 
which  was  nearly  singular,  the  approximate  solution  would  gen- 

erally be  of  large  norm  and,  quite  likely,  difficult  to  obtain  accurately. 

A small  approximate  would,  however,  result  if  we  kept  physical  coor- 

dinates known  to  be  primarily  slowly-varying  in  x^,  with  coordinates 

dominated  by  fast  boundary-layer  variations  being  kept  in  x_.  Such  a 

Y1  ' 

splitting  of  x would  correspond  roughly  to  the  splitting  y = 

y2 

into  slow  and  fast  components.  When  x^  and  are  weakly  coupled  like 

this,  the  need  for  the  decoupling  matrix  L in  (9)  is  lessened,  so  we  would 
expect  L to  be  small.  Formally,  then,  it  is  natural  to  introduce  a slow- 
mode coupling  ratio 


-15- 


(38)  p = min  (IM..I/IM.  «) 

s 1 1 II 

where  the  minimization  takes  place  over  all  possible  reorderings  of  the 
rows  of  the  matrix  M.  If  we  could  obtain  = 0,  and  x^  would 

be  decoupled  and  we'd  have  L = 0 by  (34).  The  practical  significance 
of  such  reorderings  can  be  seen  in  our  numerical  results  and  is  stressed 
in  the  aircraft  dynamics  examples  of  Teneketzis  and  Sandell  (1977),  though 
it  is  certainly  not  necessary  when  the  original  matrix  is  nonsingular. 

An  efficient  method  of  improving  the  initial  approximation  to 

the  decoupling  matrix  L can  be  obtained  by  using  a successive  approxi- 
mations scheme  on  the  algebraic  Riccati  equation  (16)  with  as  the 

first  iterate.  Using  (11),  we  can  rewrite  (16)  as 

(39)  L = b”1  (LB;l  + A21  - LA12L)  = B21(LAn  + Ay). 

Thus,  we  naturally  define  the  iterates 

(40)  Lj+1  = (A22  + LjA12)"1(LjAn  + A^),  j = 0,  1,  .... 

expecting  them  to  converge  to  the  unique  L desired.  This  iteration 
scheme  features  rapid  convergence  provided  is  close  to  L and  e is 

small  enough  (cf.  Kleinman  (1968)  who  discussed  the  analogous  Hamiltonian 
problem).  Indeed,  a proof  of  convergence  follows  from  a contraction  map- 
ping argument.  If  we  introduce 


(41) 
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we  find  that  (40)  is  equivalent  to 

(42)  D = (A22  + LjA12)"1Rj,  j - 0,  1,  .... 


for  the  residual 


(43) 


R.  = 
J 


LjAll 


' A22Lj  “ 


LjA12Lj 


+ A 


21' 


The  procedure  terminates  when  one  finds  a suitably  small  R . (In  par- 
ticular, if  Rq  is  judged  small  enough,  we  wouldn't  begin  the  iteration 
scheme  (40)  or  (42).)  Further  manipulating  (42)  and  (43)  implies  that 


(44) 


Jj+1  = (A22  + LjA12)  Dj (A11  ~ A12Lj+l) ' 


We  note  that  A 22  + L^A.^  will  be  a good  approximation  to  B7,  if 
is  a good  approximation  to  L.  Then,  the  inverse  used  in  (40)- (44)  is 
legitimate.  Likewise,  A^  - A^L  ^ then  be  a good  approximation  to 

B^,  and  since  IB  I/HB^I  = 0(e),  (44)  implies  that  our  scheme  ultimately 

has  a rate  of  convergence  of  order  e.  A closely  related  approach  is 
presented  in  Kokotovic  (1975)  and  Chow  and  Kokotovic  (1976).  Instead  of 
(37),  they  used  the  initial  iterate  * A22A21’  appropriate  if  A22  is 
nonsingular  and  if  the  system  is  sufficiently  decoupled  so  that  L is 
small  and  thereby  easy  to  compute  (cf.  (40)).  Thus,  their  convergence 
criterion  is  unnecessarily  restrictive.  Numerical  experiments  indicate 
that  our  scheme  is  successful  even  when  a rather  crude  solution  to  (37)  is 
used  for  L^,  and  when  e is  not  too  small. 

After  numerically  obtaining  the  decoupling  matrix  L,  we  need  to 
solve  the  linear  equation  (17)  for  the  matrix  K.  Rewriting  it  as 
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(45)  K - (B.jK  + A12)B2X 
suggests  the  iteration  scheme 

(46)  K = (Bf.  + A12)B21  K0  = 0,  j = 0,  1,  2,  .... 

converging  to  its  unique  solution.  A straightforward  argument  implies  a 
rate  of  convergence  of  order  e.  Here,  unlike  for  (40),  the  initial 
iterate  is  not  crucial.  Alternatively,  one  can  use  the  representations 
(32)  and  (35)  to  show  that 

(47)  K = -M12Q22  = MnQ12 

(cf.  Anderson  (1978)),  though  this  requires  more  information  about  the 
modal  decomposition  of  A than  (46). 

4.  Implementation 

We  propose  to  obtain  a reduced  order  matrix  approximation  for  two-time- 
scale  systems  x = Ax  via  the  following  numerical  algorithm: 

1)  Obtain  approximate  eigenvalues  of  A and  order  then  by  increasing 
moduli  into  sets  S and  F satisfying  (7).  Note  that  different 
choices  are  possible  for  the  number  n^  of  slow  eigenvalues  for 
any  given  matrix  A. 

2)  Determine  an  n * dimensional  real  matrix 

columns  span  the  right  eigenspace  of  A corresponding  to  the  slow 
eigenvalues.  (If  n^  > n/2,  an  alternative  based  on  left  eigen- 
vectors is  suggested  (cf.  (36)). 


M 


11 
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3)  (optional)  Reorder  the  variables  x so  that  is 

minimized. 

4)  Find  an  approximate  solution  L^  of  LgMjj  ° -M21  ^or  *ts 
alternate) . 

5)  (optional)  Improve  the  accuracy  of  the  approximation  by 

iteration.  (a)  Stop  if  II  Rj  is  less  than  a prescribed  toler- 
ance for  Rj  = L^A^  + A2l  - LjA12Lj  ~ A22Lj"  ^)  Set  Lj+1  = 

+ (A^2  + LjA^)  V and  return  to  (a)  with  j = j + 1. 

6)  The  approximate  reduced  order  fundamental  matrix  (appropriate 
within  (0,1))  is 


x(j)  _ f "l  | e(All”A12Lj)t: 


where  j is  zero  or  the  last  index  used  in  5) . 

This  algorithm  has  been  applied  to  a number  of  physical  systems  with 
order  n ranging  up  to  32.  Two  different  programs  were  used  for  the 
eigenanalysis  of  steps  1)  and  2).  They  are  the  EISPACK  subroutines,  which 
are  available  without  cost  from  the  Applied  Mathematics  Division  of  the 
Argonne  National  Laboratory  (cf.  Smith  et  al.  (1976)),  and  EIGRF,  which  is 
part  of  the  proprietary  subroutines  sold  by  the  International  Mathematics 
and  Statistics  Library  (IMSL)  of  Houston.  The  procedure  was  implemented 
on  the  64  bit  CYBER  175  and  the  36  bit  DEC  10  computers  at  the  University 
of  Arizona  Computing  Center. 


Example  1.  This  fourth  order  model  of  F-8  aircraft  longitudinal 
dynamics  (cf.  Etkin  (1972)  and  Teneketzis  and  Sandell  (1977))  has  slow 
eigenvalues  s „ - -0.0075  ± 10.076  and  fast  eigenvalues  f - - -0.94  ± 

1 9 Z L 9 l 

i3.0,  corresponding  to  the  small  parameter  e = * 0.024.  The 


L9- 


physical  variables  involved  are  velocity  variation  (ft. /sec.)*  flight 

path  angle  x2  (rad.),  angle  of  attack  x^  (rad.),  and  pitch  rate  x^ 
(rad. /sec.).  The  first  two  are  primarily  slow,  while  the  latter  two  are 


predominantly  fast.  The  four  different  orderings  of  the  x coordinates 

with  n^  = n2  = 2 show  that  the  ratio  /HM^H  and  the  corresponding 

-4  -3 

II LU  are  loosely  related.  The  numbers  are  (1.9  * 10  , 1.1  x 10  ), 

(2.4  x lo“3,  1.3  x io3),  (5.3  x 103,  6.3  x 103) , (4.2  x 102,  2.0  x 104) 

-4 

so  p =*  1.9  x 10  for  the  ordering  suggested  above.  The  success  of  the 


iterative  procedure  is  illustrated  in  Figure  1 where  log  II Rj  is  plotted 
against  j.  We  find  that  HR  /lR  I r&  t until  the  (machine-dependent) 


Figure  1. 

The  convergence  of  L for  the  F-8  aircraft  longitudinal  dynamics  model. 
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maximum  possible  numerical  accuracy  is  achieved.  These  results  also  show 


that  the  approximate  matrix 


M 


11 


M 


21 


obtained  using  IMSL  was  better  than 


with  EISPACK,  since  the  resulting  was  more  accurate.  Further  tests 

show,  however,  that  considerable  eigenvector  inaccuracies  may  be  tolerated 
with  the  algorithm  still  converging  to  the  desired  solution  when  iteration 
is  used.  Indeed,  the  algorithm  converges  whatever  ordering  is  used  for  the 
x components.  Most  important,  however,  is  the  fact  that  trajectories  for 
initial  value  problems  are  well  approximated,  outside  an  initial  boundary 
layer  region,  by  the  reduced  model.  Using  initial  conditions  (x^(0),  x^(0), 
x^(0),  x^(0))  = (100,  0,  0,  0.5),  the  aircraft  response  for  the  predomi- 
nantly slow  and  fast  variables  x^  and  x^  are  shown  in  Figures  2 and  3, 
both  for  the  full  fourth  order  model  and  the  reduced  second  order  model. 

After  the  first  few  seconds  the  responses  are  indistinguishable,  so  they  are 
not  pictured. 


1 

(ft/sec) 


Figure  2. 

Velocity  (slow  variable)  vs.  time. 


Figure  3. 

Angle-of-attack  (fast  variable)  vs.  time. 


1 

Example  2.  This  16  order  model  of  an  F-100  turbofan  engine  was  used 
as  the  theme  problem  for  the  recent  International  Forum  on  Alternatives  for 
Multivariable  Control  (cf.  Sain  (1977)).  Two  of  the  x coordinates  are 
shaft  speeds,  three  are  pressures,  and  eleven  are  temperatures.  The  eigen- 
values (ordered  in'  magnitude)  are  -0.648,  -1.906,  -2.619,  -6.715  ± 

il. 312,  -17.8  ± i4. 78,  -18.6,  -21.3  ± i0.822,  -38.7,  -47.1,  -50.7, 

-59.2,  -175.7,  and  -577.0.  If  we  take  n^  = 15,  5,  and  3,  we  get 
e values  0.304,  0.371,  and  0.383,  respectively.  The  choice  n^  = 15 

is  ruled  out  because  the  reduction  in  dimensionality  isn't  significant. 

Even  though  the  remaining  e's  aren't  very  small,  our  computational 


Iteration,  i 


Figure  4. 

The  convergence  of  for  decoupling  the  three  slowest 

modes  for  the  F-100  turbofan  engine. 


-23- 


1 


References 

1.  L.  Anderson,  "Decoupling  of  two-time-scale  linear  systems,"  Pro- 
ceedings, 1978  Joint  Automatic  Control  Conference. 

2.  R.  Bellman,  Introduction  to  Matrix  Analysis,  Second  Edition,  McGraw- 
Hill,  New  York,  1970. 

3.  A.  J.  Calise,  "Singular  perturbation  methods  for  variational  problems 
in  aircraft  flight,"  IEEE  Trans,  on  Automatic  Control  21  (1976), 
345-353. 

4.  J.  H.  Chow,  Singular  Perturbation  of  Nonlinear  Regulators  and  Systems 
with  Oscillatory  Modes,  Ph.D.  thesis,  University  of  Illinois,  Urbana- 
Champaign,  1978. 

5.  J.  H.  Chow  and  P.  V.  Kokotovic,  "Eigenvalue  placement  in  two-time- 
scale  systems,"  Proc.  IFAC  Symposium  on  Large-Scale  Systems,  Udine, 
1976,  321-326. 

6.  J.  H.  Chow  and  P.  V.  Kokotovic,  "A  decomposition  of  near-optimum 
regulators  for  systems  with  slow  and  fast  modes,"  IEEE  Trans,  on 
Automatic  Control  21  (1976),  701-705. 

7.  R.  H.  Cole,  Theory  of  Ordinary  Differential  Equations,  Appleton- 
Century-Crof ts.  New  York,  1968. 

8.  W.  A.  Coppel,  Stability  and  Asymptotic  Behavior  of  Differential  Equa- 
tions , D.  C.  Heath,  Boston,  1965. 

9.  W.  A.  Coppel,  "Dichotomies  in  stability  theory,"  Lecture  Notes  in  Math. 
629,  Springer-Verlag,  Berlin,  1978. 

10.  E.  J.  Davison,  "A  method  for  simplifying  linear  dynamic  systems,”  IEEE 
Trans.  Automatic  Control  11  (1966),  93-101. 

11.  W.  E.  Ferguson,  Jr.,  A Singularly  Perturbed  Linear  Two-Point  Boundary 
Value  Problem,"  Ph.D.  Dissertation,  California  Institute  of  Technology, 
Pasadena,  1975. 

12.  J.  E.  Flaherty  and  R.  E.  O'Malley,  Jr.,  "The  numerical  solution  of 
boundary  value  problems  for  stiff  differential  equations,"  Math,  of 
Computation  31  (1977),  66-93. 

13.  G.  H.  Golub  and  J.  H.  Wilkinson,  "Ill-conditioned  eigensystems  and  the 
computation  of  the  Jordan  canonical  form,"  SIAM  Review  18  (1976), 
578-619. 

14.  W.  A.  Harris,  Jr.,  "Singularly  perturbed  boundary  value  problems 
revisited,"  Lecture  Notes  in  Math.  312,  Springer-Verlag,  Berlin,  1973, 
54-64. 

15.  P.  W.  Hemker,  A Numerical  Study  of  Stiff  Two-point  Boundary  Problems, 
Ph.D.  Dissertation,  University  of  Amsterdam,  1977. 


-24- 


16.  M.  W.  Hirsch  and  S.  Smale,  Differential  Equations.  Dynamical  Systems, 
and  Linear  Algebra,  Academic  Press,  New  York,  1974. 

17.  F.  C.  Hoppensteadt  and  W.  L.  Miranker,  "Differential  equations  having 
rapidly  changing  solutions:  Analytic  methods  for  weakly  nonlinear 
systems,"  J.  Differential  Equations  22  (1976),  237-249. 

18.  D.  L.  Kleinman,  "On  an  iterative  technique  for  Riccati  equation 
computation,"  IEEE  Trans.  Automatic  Control  13  (1968),  114-115. 

19.  P.  V.  Kokotovic,  "A  Riccati  equation  for  block-diagonalization  of  ill- 
conditioned  systems,"  IEEE  Trans.  Automatic  Control  20  (1975),  812-814. 

20.  P.  V.  Kokotovic,  R.  E.  O'Malley,  Jr.,  and  P.  Sannuti,  "Singular  pertur- 
bations and  order  reduction  in  control  theory  - An  overview,"  Automatica 
12  (1976),  123-132. 

21.  H.-O.  Kreiss,  "Difference  methods  for  stiff  ordinary  differential 
equations,"  SIAM  J.  Num,  Anal.  15  (1978),  21-58. 

22.  V.  Kucera,  "A  review  of  the  matrix  Riccati  equation,"  Kybernetika  9 
(1973),  42-61. 

23.  W.  L.  Miranker,  The  Computational  Theory  of  Stiff  Differential  Equations, 
Instituto  per  le  Applicazioni  del  Calcolo  "Mauro  Picone,"  Rome,  1975. 

24.  C.  B.  Moler  and  C.  F.  Van  Loan,  "Nineteen  ways  to  compute  the  exponen- 
tial of  a matrix,"  SIAM  Review  20  (1978). 

25.  N.  Narasimhamurthi  and  F.  F.  Wu,  "On  the  Riccati  equation  arising  from 
the  study  of  singularly  perturbed  systems,"  Proceedings  1977  Joint 
Automatic  Control  Conference,  1244-1247. 

26.  R.  E.  O'Malley,  Jr.,  "Singular  perturbations  and  optimal  control," 
Lecture  Notes  in  Math. 

27.  M.  K.  Sain,  "The  theme  problem,"  Proceedings  International  Forum  on 
Alternatives  for  Multivariable  Control,  1977,  1-12. 

28.  C.  A.  Skira  and  R.  L.  De  Hoff,  "A  practical  approach  to  linear  model 
analysis  for  multivariable  turbine  engine  control  design,"  Proceedings 
International  Forum  on  Alternatives  for  Multivariable  Control,  1977, 
29-44. 

29.  B.  T.  Smith,  J.  M.  Boyle,  J.  J.  Dongarra,  B.  S.  Garbo,  Y.  Ikeke,  V.  C. 
Klema,  and  C.  B.  Moler,  "Matrix  eigensystem  routines-EISPACK  guide," 
Second  Edition,  Lecture  Notes  in  Computer  Science  6,  Springer-Verlag, 
Berlin,  1976. 

30.  G.  W.  Stewart,  Introduction  to  Matrix  Computations,  Academic  Press, 

New  York,  1973. 

31.  D.  Teneketzis  and  N.  R.  Sandell,  Jr.,  "Linear  regulator  design  for 
stochastic  systems  by  a multiple  time-scales  method,"  IEEE  Trans. 
Automatic  Control  22  (1977),  615-621. 


-25- 


32.  J.  E.  Van  Ness,  "Methods  of  reducing  the  order  of  power  system  models 
in  dynamic  studies,"  Proceedings,  1977  IEEE  International  Symposium  on 
Circuits  and  Systems,  858-863. 

33.  R.  Willoughby  (editor),  Stiff  Differential  Systems,  Plenum  Press,  New 
York,  1974. 

34.  F.  F.  Wu  and  N.  Narasimhamurthi,  A New  Algorithm  for  Modal  Approach  to 
Reduced  Order  Modeling,  Report  ERL-M13,  Engineering  Research  Laboratory, 
University  of  California,  Berkeley,  1976. 


Unclassified 

SECURITY  CLASSIFICATION  OF  THIS  PAGE  (Whtn  Dmtm  Entered; 

r REPORT  DOCUMENTATION  PAGE 


1.  REPORT  NUMBER  2.  GOVT  ACCESSION  NO. 

Twenty-One 


4.  TITLE  (and  Subtitle) 

Singular  Perturbations,  Order  Reduction, 
and  Decoupling  of  Large  Scale  Systems 


7.  authorcs; 


R.  E.  O'Malley,  Jr.  and  L.  R.  Anderson 


9 PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

Department  of  Mathematics  ' 
University  of  Arizona 
Tucson,  Arizona  85721 


111.  CONTROLLING  OFFICE  NAME  ANO  ADDRESS 

Mathematics  Branch 

I Off ice  of  Naval  Research 

Arlington,  Virginia  


«.  MONITORING  AGENCY  NAME  8 ADDRESSf//  dlttorent  from  Controlling  Olllce) 


16.  DISTRIBUTION  STATEMENT  (ot  this  Report) 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


3.  RECIPIENT'S  CATALOG  NUMBER 


5.  TYPE  OF  REPORT  6 PERIOD  COVERED 

Technical  Report 
May  1978 


6.  PERFORMING  ORG.  REPORT  NUMBER 


8.  CONTRACT  OR  GRANT  NUMBERf'e; 


N00014-7  6-C-0326 


10.  PROGRAM  ELEMENT.  PROJECT,  TASK 
AREA  A WORK  UNIT  NUMBERS 


NR041-466 


12.  REPORT  OATE 

May  1978  


13.  NUMBER  OF  PAGES 

Twenty-Five  (25) 


IS.  SECURITY  CLASS.  ( of  thlt  report) 

Unclassified 


15«.  DECLASSI  FI  CATION/ DOWNGRADING 
SCHEDULE 


Approved  for  public  release.  Distribution  unlimited. 


19.  KEY  WOROS  (Continue  on  reverse  side  II  necessary  and  Identity  by  block  number) 


Singular  perturbations,  stiff  equations,  reduced-order  modeling, 
modal  decoupling,  large  scale  systems 


2o\l ABSTRACT  (Continue  on  reverse  side  it  necessary  and  Identify  by  block  number) 


Singular  perturbation  methods  have  recently  provided  new  techniques  for 
optimal  control  and  stability  problems.  Their  primary  advantage  is  the 
order  reduction  achieved  when  small  parameters  are  set  equal  to  zero.  We 
present  an  approach  for  identifying  small  parameters  by  separating  slow 
and  fast  modes.  A computational  procedure  is  developed  for  simple  linear 
systems  and  tested  numerically.  Connections  to  current  engineering 
literature  is  stressed. 


DD  ,5S“»  1473  EDITION  OF  I NOM  68VS  OBSOLETE 


S/N  0102  LF  014  6601 


CD  - V 


Unclassified 

SECURITY  CLASSIFICATION  OF  THIS  PAGE  (Whon  Dmtm  tntmrmd) 


